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Abstract 

These  lectures  review  some  basic  results  of  detonation  theory  and 
are  specialized  to  serve  as  an  introduction  to  the  theory  of  Detona¬ 
tion  Shock  Dynamics  (DSD).  The  theory  of  DSD  is  a  time-dependent, 
multi-dimensional  theory  for  the  propagation  of  near-Chapman-  Jouguet 
(CJ)  detonations.  This  theory  is  especially  relevant  to  applications  of 
detonation  propagation  in  condensed  explosives  and  the  basic  dynam¬ 
ics  in  freely  propagating  gaseous  explosives  as  well.  The  theory  that 
we  present  is  based  on  rigorous  mathematical  arguments  and  rational 
approximations  for  an  assumed  model  of  the  explosive  material.  The 
material  is  described  by  the  compressible,  reactive  Euler  equations 
for  a  given  equation  of  state  and  kinetic  rate  law  for  the  release  of 
exothermic,  chemical  energy.  This  theory  is  the  basis  for  the  engi¬ 
neering  ” Method  of  Detonation  Shock  Dynamics”,  that  is  now  being 
used  in  the  design  of  explosive  systems.  Lectures  on  the  Method  of 
Detonation  Shock  Dynamics  are  being  planned  as  a  sequel.  This  first 
set  of  lectures  were  given  in  the  summer  of  1992  at  Eglin  Air  Force 
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1  Outline  of  the  Lectures 

The  outline  of  theses  lectures  is  eis  follows: 

•  A  brief  description  of  Detonation  Shock  Dynamics 

•  The  basic  equations 

-  The  compressible,  reacting  Euler  equations 

—  The  Rankine-Hugoniot  (RH)  relations  in  shock  normal  coordi¬ 
nates 

—  Intrinsic  shock  dynamics  of  the  RH  relations  represented  in  Carte¬ 
sian  coordinates 

•  Review  of  steady,  ID,  Zeldovich,  Neumann,  Doering  (ZND)  theory 

•  Shock  attached  coordinates 

-  General  form  of  the  governing  equations 

•  The  detonation  shock  speed  (Dn)  -  curvature  (k)  relationship 

—  Derivation  of  the  "master”  equation  on  a  central  streamline  nor¬ 
mal  to  the  shock  surface 

-  Derivation  of  the  D„  —  k  relationship  under  the  assumptions  of 
small  shock  curvature,  (/c  <<  1),  and  slow  time  variations  {djdt  « 
1) 

In  the  sequel  to  these  lectures  we  plan  to  add  the  following  topics; 

•  Remarks  on  conditions  near  interfaces  and  suitable  boundary  condi¬ 
tions 

•  Properties  and  solution  of  the  Dn  —  n  relations.  The  "Method  of  DSD” 

—  Basic  dynamics  of  a  DSD  wave 

-  Comparison  with  numerical  experiment 

-  Comparison  with  physical  experiment 
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•  Research  issues 

—  Extension  of  the  theory  at  interfaces  and  for  wave  -  wave  interac¬ 
tions 

—  Hydrodynamic  stability 
—  General  equation  of  state 

•  Constraint  of  the  rate  equation  by  matching  with  the  Dn  ~  ^ 
relation 

•  Imbedding  of  DSD  waves  in  ’’fast  -  time”  dynamics 

2  Introduction:  A  brief  description  of  Deto¬ 
nation  Shock  Dynamics 

Detonation  Shock  Dynamics  (or  DSD)  is  the  name  given  to  a  body  of  multi¬ 
dimensional  theory  that  describes  the  dynamics  of  ”  near-  Chapman- Jouguet” 
detonations  and  is  named  after  Whitham’s  theory  of  "Geometrical  Shock 
Dynamics”  because  of  the  similarity  of  the  mathematical  structure  of  the 
theories. 

The  basic  result  of  DSD  theory  is  that  under  suitable  conditions,  the 
detonation  shock  in  the  explosive  propagates  according  to  the  simple  formula 

Dn  =  £>CJ  -«(«),  (1) 

where  Dn  is  the  normal  velocity  of  the  shock  surface,  Dcj  is  the  ID,  steady, 
Chapman-Jouguet  velocity  for  the  explosive,  and  a(K)  is  a  function  of  cur¬ 
vature  K,  that  is  determined  by  the  material  properties  of  the  explosive.  A 
sketch  of  a  typical  Dn  —  n  relation  is  shown  in  figure  1 . 

The  curvature  used  in  this  formula  is  the  local  shock  curvature  at  given 
point  on  the  shock,  with  a  corresponding  attached  normal.  The  outward 
normal  direction  is  defined  to  be  positive  when  it  is  pointing  in  the  direction 
of  the  unreacted  explosive.  Therefore  when  the  curvature  is  positive,  k  >  0 
the  shock  shape  is  convex  and  the  detonation  is  expanding  and  its  surface 
area  is  increasing.  Negative  curvature,  «:  <  0,  corresponds  to  a  converging 
detonation  with  a  concave  shape.  Figure  2  shows  these  cases  for  typical 
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Figure  1:  The  Dn  —  k  relation  for  a  typical  condensed  phase  explosive  after 
Bdzil  et.  al.’s  calibration  of  PBX  9502,  [1] 

detonation  shocks.  The  detonation  shock  samples  the  curvature  of  the 
relation  and  representative  points  on  the  detonation  shock  are  indicated  at 
time  t  =  0  (say)  by  a  broken  curve.  The  solid  line  shows  the  advancement  of 
the  shock  at  time  t  +  At. 

2.1  Surface  evolution  defined  by  the  Dn  —  k  relation 

The  £)„  —  K  relation  is  a  surface  evolution  equation  for  the  detonation  shock 
surface.  Once  a  coordinate  system  is  chosen,  the  —  «  relation  can  be 
shown  to  be  a  nonlinear,  parabolic  partial  differential  equation.  A  simple 
but  low  resolution  solution  for  the  motion  of  the  surface  can  be  obtained 
graphically.  One  would  discretize  the  surface,  marking  many  points,  at  equal 
arclength  (say)  on  the  shock  at  time  t,  draw  the  local  normals  to  the  surface, 
measure  the  sum  of  the  principle  radius  of  curvatures  ,  k  =  {k\  +  K2)  each 
point  on  the  shock  surface.  Then  one  would  look  up  the  appropriate  value  of 
Dn  for  the  measured  curvature  k,-,  and  then  mark  a  new  point  on  the  normal 
at  a  distance 


An,-  =  D„(/c,)At.  (2) 

The  advanced  surface  then  is  approximated  as  that  one  that  passes  through 
all  of  the  new  points  advanced  An,-  along  their  respective  shock  normals.  This 
simple  algorithm  is  equivalent  to  a  discretized,  first-order  accurate  numerical 
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Diverging  Geometry 


Converging  Geometry 


Figure  2:  Figure  2a.  shows  a  diverging  portion  of  a  near-CJ  detonation  shock. 
Figure  2b.  shows  a  converging  portion  of  a  near-CJ  detonation  shock.  The 
values  of  the  curvature  at  the  numbered  points  are  indicated  on  the  Dn  —  k 
curve  shown  in  figure  1. 

solution  for  the  surface  evolution  equation,  which  is  a  second  order  partial 
differential  equation  in  space  and  first  order  in  time.  Figure  2  indicates  the 
graphical  surfece  construction  for  a  converging  and  diverging  cases.  The 
broken  curve  indicates  the  shock  surface  at  the  current  time,  t,  and  the  solid 
curve  indicates  the  shock  surface  at  a  previous  time  t  +  At. 

Since  the  Dn-  k  relationship  is  usually  monotonically  decreasing  with 
increasing  curvature,  regions  of  the  detonation  shock  that  are  diverging  are 
slowed  and  regions  of  the  shock  that  are  converging  are  sped  up,  relative  to 
Dcj- 


2.2  Importance  of  the  Dn  —  «  relationship  for  multi¬ 
dimensional  theory  and  application 

Here  we  mention  a  few  important  points  about  the  properties  of  the 
relationship  and  their  implications  for  basic  detonation  theory  and  explosives 
engineering. 
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•  The  Dn  —  n  curve,  combined  with  an  analysis  of  the  explosive  con¬ 
finement  provided  by  the  walls  of  the  stick,  predicts  the  well-known 
diameter  effect  in  rate  stick  geometry.  (The  diameter  effect  reflects  the 
fact  that  the  steady  axial  velocity  of  curved  detonation  is  a  function 
that  decreases  as  the  stick  radius  decreases.) 

•  DSD  theory  shows  that  the  —  k  relationship  is  an  intrinsic  (i.e. 
material)  property  of  the  explosive.  Therefore  the  Dn  —  k  curve  can 
be  measured  in  one  geometry  and  then  be  used  to  make  predictions  in 
other  geometries. 

•  The  Dn  —  n  relationship  is  a  correction  of  the  Huygens  construction 
rule  that  has  traditionally  been  used  to  describe  ideal  detonation  prop¬ 
agation.  In  the  Huygens  construction,  the  detonation  normal  velocity 
is  constant  and  equal  to  its  CJ  value.  The  Z)„  —  k  relationship  is  found 
theoretically  to  depend  on  reaction  zone  chemistry  and  describes  the 
modification  of  the  normal  detonation  shock  velocity  from  CJ.  For  ex¬ 
ample,  the  longer  the  ID  reaction  zone  is,  the  more  important  the 
relative  reaction  zone  effect  is  on  the  dynamics  of  that  explosive.  This 
is  an  important  consideration  for  insensitive  high  explosives  that  typ¬ 
ically  have  long  reaction  zones  compared  to  sensitive  high  explosives 
with  short  reaction  zones. 

•  The  shock  surface  evolution  can  be  carried  out  with  a  small  computa¬ 
tional  effort,  when  it  is  applicable,  relative  to  that  of  direct  numerical 
simulation.  Prediction  of  the  motion  of  the  detonation  shock  surface  by 
DSD,  requires  the  solution  of  a  parabolic  PDE.  The  same  simulation 
carried  out  using  a  Direct  Numerical  Simulation  (DNS)  of  the  Euler 
equations,  requires  a  solutipn  of  a  minimum  of  five  (in  2D),  or  six  (in 
3D)  coupled,  hyperbolic  equations.  A  more  important  consideration  is 
the  that  the  size  of  the  simulation  done  for  DNS  is  severely  restricted 
because  of  accuracy  limitations.  DSD  calculations  are  not  limited  in 
the  same  way  and  can  be  used  on  much  larger,  device-sized  regions, 
without  loss  of  accuracy,  within  the  assumptions  of  the  theory. 

•  The  Dn  —  n  relationship  is  sensitive  to  the  form  of  the  rate  law  (chem¬ 
istry)  and  therefore  may  be  an  excellent  experimental  indicator  and  a 
way  to  measure  the  average  reaction  zone  chemistry. 
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2.3  A  Brief,  incomplete  history  of  the  Dn  —  k.  relation 

Here  we  outline  some  of  the  most  important  contributions  to  the  development 
of  the  theory.  This  list  is  incomplete  but  indicates  the  main  theoretical 
contributions. 

The  history  of  the  Dn  —  /c. relation  starts  with  Eyring  et.  al.  (1949),  [2] 
who  suggested  that  the  curvature  of  the  detonation  wave  might  have  a  sig¬ 
nificant  effect  on  its  propagation.  The  next  important  contribution  was  that 
of  Wood  and  Kirkwood  (1954),  [3]  who  attempted  to  derive  the  relationship 
on  a  central  streamtube.  Their  work  was  based  on  ad  hoc  assumptions,  how¬ 
ever  they  outlined  the  central  analysis.  Fickett  and  Davis  generalized  the 
discussion  of  the  problem  posed  by  Wood  and  Kirkwood  [14]  in  a  section  of 
their  book  entitled  ’’Slightly  Divergent  Flow”  and  used  the  relative  velocity 
squared  and  the  reaction  progress  variable,  phase-plane.  Bdzil  (1981),  [4] 
made  a  critical  contribution  that  is  the  origin  of  the  asymptotic  theory  of 
DSD.  In  his  paper,  Bdzil  derived  the  diameter  effect  that  relates  the  axial 
speed  of  the  detonation  to  the  diameter  of  the  rate.  To  do  this,  he  used  formal 
asymptotic  arguments  based  on  the  small  reaction  zone  thickness  compared 
against  the  diameter  of  the  stick  and  established  an  ordinary  differential 
equation  for  the  steady  shock  shape  across  the  stick.  The  steady  Dn  —  k 
relationship  is  implicit  in  Bdzil’s  paper.  The  first  asymptotic  derivation  of 
the  3D,  unsteady  D„  —  k  relation,  was  derived  (and  presented  essentially  in 
the  form  discussed  in  these  notes)  by  Stewart  and  Bdzil  (1988),  [5]. 

Other  contributions  to  the  asymptotic  theory,  experimental  and  numer¬ 
ical  implementation  of  DSD,  include  those  of  Lambourn,  [6],  Bukiet  and 
Menikoff,  [7],  and  Klein  and  Stewart,  [8].  There  have  been  other  contribu¬ 
tions  that  use  ad  hoc  approximations  that  we  don’t  mention  here.  Bdzil  has 
pioneered  the  practical  application  of  DSD  for  engineering  applications,  [9] 
and  Lambourn  and  Swift  have  worked  on  the  engineering  applications  in  a 
similar  effort  in  England  [10].  A  list  of  related  references  is  found  in  the 
bibliography. 


3  Governing  Equations  and  the  Basic  Model 

The  basic  mathematical  model  of  explosive  materials  is  the  compressible  Eu¬ 
ler  equations  with  reaction.  A  state  variable  theory  is  assumed,  where  the 
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explosive  is  described  by  its  thermomechanical  properties  for  the  bulk  ma¬ 
terial.  The  basic  mechanical  variables  are  the  velocity,  u,  the  density  p  and 
the  thermodynamics  pressure  p.  Chemistry  is  modeled  directly  in  the  equa¬ 
tion  of  state  by  introducing  additional  thermodynamic  degree(s)  of  freedom 
that  measures  the  reaction  progress  of  a  net,  forward,  exothermic  chemical 
reaction.  One  global  progress  variable.  A,  can  be  used,  but  a  reaction  scheme 
can  also  be  considered  where  A  is  replaced  by  a  vector  A,-.  Specification  of 
an  equation  of  state  of  the  form  e(p,  p,  A)  equation  of  state,  rate  law  for  A 
(or  rate  laws  for  A,),  is  assumed  to  be  given  at  the  outset  to  describe  the 
explosive. 

3.1  Material  description 

3.1.1  Equation  of  state 

In  these  notes,  we  will  illustrate  the  DSD  theory  mainly  for  the  polytropic 
equation  of  state, 

e  =  ^^-<3A.  (3) 

/>7-  1 

Here  7  is  the  polytropic  exponent  and  Q  is  the  heat  of  combustion.  This 
equation  of  state  is  the  appropriate  one  for  a  description  of  a  gaseous  explo¬ 
sive.  The  polytropic  equation  of  state  is  often  used  to  describe  the  expansion 
of  explosive  products  by  allowing  7  to  have  an  artificially  higher  value  than 
that  usually  allowed  for  gases,  i.e.  7  ~  2.5  —  3.  This  EOS  also  has  the  ad¬ 
vantage  that  a  relatively  large  body  of  theoretical  results  exist  for  it.  These 
include  resolved  numerical  and  parametric  studies  of  the  DSD  theory  and 
hydrodynamic  stability,  [11],  [12].  DSD  theory,  however  is  not  restricted  to 
ideal  EOS  and  nonideal  explosive  equations  of  state  can  be  used. 

For  example,  one  might  consider  a  typical  model  of  an  explosive  that  in¬ 
terpolates  between  an  equation  of  state  that  matches  the  near-shock  EOS, 
when  A  ~  0  and  matches  to  the  CJ  products  EOS  when  A  ~  1  (near  complete 
reaction).  The  near-shock  EOS,  typically  of  the  Mie-Gruniesen  form,  is  con¬ 
structed  from  a  Taylor  series  expansion  in  the  vicinity  of  an  experimentally 
determined  shock  Hugoniot  for  the  unreacted  explosive.  This  form  is  then 
guaranteed  to  at  least  fit  the  unreacted  explosive’s  response  to  a  shock.  Let 
the  near-shock  EOS  be  e  =  €.he{PiP)- 
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An  example  of  the  near-complete  reaction  EOS  is  the  JWL-EOS,  which  is 
also  of  Mie-Gruniesen  form,  constructed  from  a  Taylor  series  expansion  that 
uses  the  CJ  isentrope  as  a  reference  curve,  [14].  Suppose  we  refer  to  that 
EOS  as  e  =  epR{p,p).  Using  the  reaction  progress  variable  to  provide  an 
interpolation  between  the  two  behaviors,  the  composite  EOS  can  be  written 
in  a  form  that  resembles  that  for  mixtures  of  gases,  for  example 

e  =  (1  -  X)eHE{p,  p)  +  AepR(p,  p).  (4) 

3.1.2  Burn  rate  law 

The  burn  rate  law,  is  an  evolution  equation  for  the  material  rate  of  change 
of  the  progress  variable  A  within  the  explosive  material  particle.  The  form 
of  this  rate  law,  derived  from  considerations  of  the  microstructure,  has  been 
an  uncertain  element  in  the  modeling  of  condensed  explosives.  One  clear 
constraint  is  that  the  reaction  should  stop  when  A  =  1,  which  suggests  that 
the  reaction  rate  should  be  a  function  of,  and  vanishes  with  1 —A.  The  rate  law 
is  usually  constrained  by  comparison  with  macroscopic  physical  experiment. 
A  typical  form  is 


=  r(p, p,  A)  =  fc(p,  p)(l  -  A)",  (5) 

where  DXjDt  is  the  material  derivative  of  A.  Common  forms  for  the  rate 
premultipUer  fc(p,  p)  include: 

k{p^  p)  =  fc,  a  constant,  (6) 

jk(p,p)  =  fc  ca:p(-f;/(p/p)),  an  Arrhenius  form,  (7) 

fc(p,  p)  =  k  exp{A  +  Bp  +  Cp^),  a  Forest  Fire  form.  (8) 


In  the  above,  k,  E,  A,  B  and  C  are  constants. 

In  summary,  the  explosive’s  material  properties  are  described  entirely  by 
the  e(p,  p.  A)  equation  of  state  and  the  functional  form  for  the  rate  law,  i.e., 
r(p,p,A). 
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3.2  The  conservation  laws 


The  equations  of  motion  are  the  Euler  equations  for  a  compressible,  reacting 
material.  The  material  or  Lagrangian  derivative  is  represented  as 


_D 

Dt 


^  +  V. 

dt 


(9) 


The  statements  of  conservation  of  mass,  momentum,  energy  and  statement 
and  the  statement  that  the  material  time  rate  of  change  of  the  reaction 
progress  variable  is  equal  to  the  local  reaction  rate  are  given  by 


where  v  =  If p  and 


(10) 

+  ^p  =  0, 

(11) 

De  Dv 

(12) 

II 

(13) 

These  equations  hold  everywhere  in  smooth  part  of  the  flow,  i.e.  with  smooth 
variations  of  the  field  variables  and  their  first  spatial  derivatives.  However, 
the  flow  field  is  generally  comprised  of  smooth  variations  which  are  inter¬ 
sected  by  gasdynamic  discontinuities  such  as  shocks,  contact  discontinuities 
and  rarefaction  fans.  The  detonation  wave,  in  particular  has  a  lead  shock, 
followed  by  a  smoothly  varying  reaction  zone. 


3.3  Rankine>Hugoniot  relations 

The  ZND  structure  of  a  detonation  wave  has  a  shock  wave  preceding  a  reac¬ 
tion  zone.  A  sketch  is  shown  of  the  pressure  and  progress  variable  is  shown 
in  figure  3.  The  lead  shock  wave  is  a  moving  boundary  whose  motion  and 
location  is  coupled  to  the  smooth  variations  behind  it.  If  the  upstream,  unre¬ 
acted  quiescent  state  is  known,  and  the  motion  of  the  shock  is  approximated, 
then  in  the  solution  scheme  the  shock  acts  as  a  boundary,  and  the  values  of 
the  variables  immediately  behind  the  shock  are  boundary  values. 
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Figure  3:  A  sketch  of  the  structure  of  a  ZND  detonation,  which  is  shock  wave 
followed  by  a  reaction  zone 

3.3.1  Shock  normal  coordinates 

The  dynamics  of  the  shock  and  the  mathematical  representation  of  the  sur¬ 
face,  in  different  coordinates  systems,  play  a  central  role  in  the  theory  of  DSD. 
The  starting  point  of  the  shock  surfax:e  dynamics  is  the  Rankine-Hugoniot 
shock  relations  which  conserve  mass,  momentum  and  energy  of  material  that 
flows  through  the  shock  surface. 

The  shock  is  assumed  to  be  a  smooth  surface  of  discontinuity,  with  an 
outWard  normal  to  the  surface,  h  in  the  direction  of  the  unreacted  material, 
and  with  two  independent  vectors,  in  the  surface  of  the  shock,  parallel 
to  the  tangent  plane  defined  by  the  normal.  A  convenient  choice  defines  these 
unit  vectors,  in  the  directions  of  the  principle  radius  of  curvature  of  the 
surface.  Figure  4  shows  a  2D  shock  relative  to  laboratory  coordinates.  We 
discuss  intrinsic  coordinates  at  length  in  section  4. 

We  use  a  (0)  subscript  to  denote  the  state  in  the  unreacted  material  and 
use  a  (s)  in  the  shocked  material.  Suppose  that  the  shock  surface  velocity  is 
given  at  each  point  by  the  vector  D.  Let  Un  =  u  •  fi,  Uf,  =  u  •  U,  Dn  =  D  •  n, 
define  the  normal  and  tangential  particle  velocities  and  the  normal  shock 
velocity.  The  normal  mass  flux  m  is  given  by  m  =  p{u  —  D)  •  h.  The 
momentum  flux  normal  and  tangential  to  the  shock  are  m(u  —  D)  •  h  and 
m{u—D)'i  respectively.  Then  the  standard  form  of  the  normal  shock  relation 
are  given  by. 
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Figure  4:  The  geometry  for  a  2D  detonation  shock  moving  in  a  fixed  labora¬ 
tory  coordinate  system 


[p{u  -  D)  •  n]o  =  [/)(u  -  D)  •  n],  =  m,  (14) 

\p  +  m{u-D)’h]o  =  \p  +  m{u  —  D)‘h]s,  (15) 

[m{D  —  u)  •  i\o  =  [m{D  -  u)  ‘  t\s,  (16) 

[e  +  p/p  +  l/2(u„  -  Dn)%  =  [e  +  pjp  +  l/2(«„  -  D„)^]„  (17) 

A,  =  0.  (18) 

Another  equivalent  (and  convenient)  form  for  the  shock  relations  is 

[p{Un  -  Dn)]o  =  [p{Un  -  Dr,)]s,  (19) 

p,  -  Po  =  m^vo  -  Vs),  (20) 
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for  i  =  1,  2  , 


(21) 


e,  -  eo  =  l/2(p,  +  Po){vo  -  v,),  (22) 

A,  =  0.  (23) 

For  a  given  shock  equation  of  state,  e(p,  p,  0)  and  for  a  given  D„,  one 
can  solve  explicitly  for  the  shock  state,  pa{po,P(ii  Dn),  etc.  The  solution  for 
the  shock  state  pi,(«n)»,(«fi)«>P*»'^*  of  Dn,  serves  as  the  boundary 

conditions  for  the  solution  of  the  Euler  equations.  However,  the  velocity  of 
the  shock,  Dn  is  ultimately  determined  by  compatibility  with  the  flow  behind 
the  shock  and  thus  is  part  of  the  complete  solution. 

3.3.2  Exercise:  Rankine  Hugoniot  algebra  for  the  ideal  EOS 
For  the  ideal  equation  of  state  with  reaction, 

e  =  -^^-gA,  (24) 

7  -  1  p 

define  C/„  =  Un  —  Dn  (the  relative  normal  velocity)  and  take  uq  =  0  so 
that  m  =  -poD,  from  the  Rankine  -  Hugoniot  relations,  eliminate  p  and  U 
in  favor  of  p,  and  obtain  a  quadratic  equation  for  p*  show  that  pa  has  the 
solution 


Ps  = 


7  +  1 


m‘ 
Po  +  — 
Po 


( 


1  m‘ 


:Po  - 


7+1  7  +  1  Po 


)■- 


7-1 
7  +  1 


2Tn^QX 


1/2 


.  (25) 


The  inert  shock  solution  is  found  by  setting  QX  =  0.  Show  that  when 
QA  =  0  the  two  solutions  are  the  undisturbed  state,  with  p*/po  =  1  (the 
minus  branch)  and  the  shock  state  (the  plus  branch)  with 

^  =  1  +  -  1),  (26) 
Po  7  +  1 

where 
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Figure  5:  The  Hugoniot  and  the  Rankine  line  plotted  in  the  p—v  plane.  The 
intersection  defines  the  shock  state. 


Ml  = 


Co  ’ 


where,  =  JPol  pQ. 


is  the  upstream  Mach  number  of  the  shock. 

Once  the  shock  pressure  ps  is  determined,  then  specific  volume  u,  and  the 
relative  normal  particle  velocity,  Us  are  calculated  by  successive  evaluations, 


Pa  -PO 

Vs  =Vq-  , 

(28) 

(/.  =  =. 

(29) 

The  solution  to  the  shock  relations  can  also  be  obtained  by  graphical 
solution  in  the  p  —  u  plane.  One  plots  the  Rayleigh  line,  (20)  (a  linear 
relationship  between  pressure  and  specific  volume)  and  the  Hugoniot  curve, 
(22).  The  intersection  of  Rayleigh  line  and  the  Hugoniot,  for  Ps/po  >  1 
defines  the  solution.  See  figure  5. 


3.3.3  Exercise:  Strong  shock  approximation 

The  strong  shock  approximation  is  the  limit  when  the  ratio  of  the  shock 
pressure  to  the  ambient  pressure  is  very  large,  i.e. 
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^  »  1.  (30) 

Po 

Show  that  in  this  limit  the  shock  relations  for  the  ideal  equation  of  state 
reduce  to, 


7  +  1 


Ps  = 


7  +  1 


Po^l, 


Ufi  —  ILfi  Dfi  — 


7  -  1 
7  +  1 


(31) 

(32) 

(33) 


3.4  Cartesian  coordinates 

Consider  a  lab-frame  Cartesian  coordinate  system,  (similar  to  that  shown  in 
figure  4),  where  z  is  in  the  main  (axial)  direction  of  propagation  of  the  deto¬ 
nation  shock  and  the  x  and  y  axes  are  in  the  transverse  directions.  Suppose 
the  shock  surface  is  described  by  the  surface  equation 


x{)=:  z- Dt-Za{x,y,t)  =  0,  (34) 

V  ^ 

where  D  is  a  constant  and  describes  the  nominally  steady  1-D  detonation 
velocity  {D  =  Dcj  ,  say).  Then  the  function  Zs(x,  t/,  t)  describes  the  deviation 
of  the  shock  locus  from  the  plane,  2  =  Dt.  These  coordinates  are  useful 
for  describing  detonation  in  a  geometry  typical  of  a  rate  stick,  where  the 
detonation  travels  down  the  axis  of  the  stick. 


3.4.1  Exercise:  Surface  Kinematics 

Given  the  moving  surface  as  described  above,  show  that  the  normal  n  is  given 
by 


n  =  —z — ^  = 

IVV-I 


ax  dy  ^  ^ 


[1  +  (|s?V  + 


1 1/2' 


^  Uv ;  J 

and  that  the  normal  velocity  of  the  shock  surface  can  be  written  as 


(35) 
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or 


Next,  apply  these  specific  representations  to  the  vector  form  of  the  Rankine- 
Hugoniot  relations,  with  u  =  Uxtx  +  UyCy  +  UzCz-  In  particular  you  will  need 
to  use  the  mass  conservation  condition  (which  involves  u„)  and  the  tangential 
momentum  conservation  condition  (which  involve  the  components  of  Uf,)  to 
obtain  conditions  for  the  jumps  of  Uj,  Ux  and  Uy. 

Show  that  in  the  strong  shock  limit,  for  the  polytropic  equation  of  state, 
in  Cartesian  coordinates,  the  Rankine-Hugoniot  conditions  can  he  solved  for 
the  shock  state  (s),  in  terms  of  the  function  z,(x,y,t)  that  describes  the 
motion  of  the  shock  as 


A,  =  0.  (43) 
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3.5  Dynamics  of  the  shock  interface  for  small  shock 
slope 

From  the  shock  relations,  we  can  make  observations  about  the  dynamics  of 
the  shock  when  the  shock  slope  is  small  relative  to  the  plane  of  the  axial 
direction.  To  make  rational  approximations  of  the  dynamics,  we  need  to 
select  the  axial,  reference  normal  that  defines  the  reference  tangent  plane 
and  then  discuss  the  motion  of  the  shock  relative  to  it. 

The  smallness  of  the  shock  slope  relative  to  the  reference  plane  is  given 
by  the  assumption  that 

—  0(c),  where  e  <<  1.  (44) 

ox 

The  distortion  of  the  shock  from  planar  shows  up  in  the  appearance  of  terms, 
like  dzsidx  in  the  shock  jump  relations.  We  see  that  shape  change  effects 
the  shock  state  by  an  ordered  amount.  In  particular,  the  shock  pressure  is 
changed  by  O(e^)  from  its  steady  ID  value  by  shape  effect  terms.  The  same 
is  true  for  the  particle  velocity  in  the  axial  direction,u,.  However,  the  leading 
order  correction  to  the  transverse  velocity  u*  is  0(c). 

The  shock  relations  of  the  last  section  also  show  that  unsteady  effects  of 
motion  of  the  shock  surface  are  of  the  same  order  of  magnitude  as  the  shock 
shape  effects  if  the  dynamics  of  the  shock  'motion  are  sufficiently  slow.  The 
explicit  assumption  that  shock  shape  and  shock  motion  effects  are  compara¬ 
ble  in  their  influence  on  the  shock  state  is  written  as 


(45) 

With  the  assumption  of  the  scalings  given  by  equations  (44)  and  (45), 
the  expansions  for  the  the  Rankine-Hugoniot  shock  states  for  Uj,  Ux  and  p, 
take  the  approximate  form. 


7  +  1 

(«r) 


D\l  + 


1  dz, 
D  dt 


.  dx 


2 


7+1  dx 


(46) 
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(47) 


(48) 


P< 


2 

7  +  1 


PoD'^ 


1  dzs 

D  dt 


The  expression  for  (uj,),  is  similar  to  (47)  and  p,  and  A,  are  unchanged  from 
(42)  and  (43)- 

The  conclusions  of  this  section  can  be  summarized  as  follows.  Small 
spatial  defects  on  the  shock  are  consistent  with  slowly  evolving  dynamics. 
Such  a  theory  of  approximation,  leads  to  a  quasi-  ID,  quasi-steady  description 
of  the  detonation  structure.  Small  shock  slope  corrections,  dzg/dx  ~  0(c) 
and  dynamic  shape  changes  {dzaldt)ID  O(e^)  lead  to 


•  O(e^)  changes  in  the  shock  pressure 

•  O(e^)  changes  in  the  normal  particle  velocity  along  axial  direction  ,  (e.g 

Uz) 

•  0(e)  changes  in  the  transverse  particle  velocities  (e.g.  Ux  ) 

For  this  special  equation  of  state  in  the  strong  shock  limit,  there  are  no  den¬ 
sity  changes  at  the  shock,  however  the  density  expansion  follows  the  velocity 
expansion  in  the  normal  direction  and  the  correction  to  the  density  at  the 
shock  is  O(c^).  The  ordering  scheme  directly  available  from  the  shock  re¬ 
lations  suggests  a  specific  set  of  scaling  relations  for  the  construction  of  an 
asymptotic  theory.  We  use  these  scalings  later  when  we  derive  the  master 
equation  in  section  5. 


4  Review  of  ID,  ZND  Detonation  Theory 

The  ZND  theory  of  detonation  structure  assumes  that  the  detonation  is  a 
shock  wave,  with  a  reaction  zone  that  follows  behind  the  shock  and  provides 
the  energy  that  supports  the  detonation.  From  the  governing  equations,  it  is 
possible  to  find  a  supersonic,  steady  traveling  wave  that  fits  this  description. 
The  analysis  of  the  structure  of  the  steady  ZND  is  the  first  step  in  a  more 
complete,  multi-dimensional  theory. 

To  obtain  the  equations  that  describe  the  smoothly  varying  part  of  the 
detonation  structure,  i.e.  the  reaction  zone,  one  assumes  that  the  wave  is 
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traveling  in  the  ^-direction  (say)  with  a  fixed  speed  D.  The  particle  velocities 
in  the  transverse  direction  are  assumed  to  be  zero  and  only  variations  in  z 
are  allowed.  The  wave  is  assumed  to  be  steady  in  a  traveling  wave  coordinate 
defined  by  n  =  2 —  Dt. 


4.1  Exercise:  Steady  wave  analysis 

Make  the  following  assumptions: 

u  =  u»e,,  (49) 

and  let  U  =  Uz{n)  —  D.  Substitute  these  forms  into  the  governing  equations 
and  show  that  the  following  ODE’s  result 


rrdp  dU  „ 

rrdU  dp  ^ 

de  dv  ^  ,  , 

—  +  p-r  =  0»  with  V  =  1/p, 
an  an 

dX 

which  can  be  written  in  conservative  form  as 


(50) 

(51) 

(52) 

(53) 


^(pU)  =  0,  (54) 

^{p  +  pU^)=0,  (55) 


— (e  +  pI  P  +  ll2U^)  —  0, 


(56) 


(57) 


(The  last  equation  of  course,  is  not  a  conservation  statement,  but  would  be 
if  there  were  no  reaction.) 
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(58) 


Integration  of  three  of  the  last  four  equations  obtains 

pU  =  Cr, 

p  +  pU^  =  C,,  (59) 

e(p,^,A)  +  ^  +  V  =  C3,  (60) 

P  2 

which  are  valid  for  all  n  within  the  reaction  zone.  In  particular  these  equa¬ 
tions  must  hold  at  the  shock,  n  =  0.  Hence  by  evaluating  these  relations  at 
the  shock,  and  in  turn,  by  using  the  shock  relations  themselves,  one  finds  (in 
the  strong  shock  approximation), 

C,  =  pM,  =  -PoD,  (61) 

Ci  =  p,  +  i>,U^  =  PoD‘,  (62) 

C3  =  e(p„/>„0)  +  l/2l/f +p./^..  =  l/2D».  (63) 

The  task  is  to  solve  these  relations  for  p,p  and  U  as  functions  of  A  and 
D,  which  in  turn  hold  throughout,  the  reaction  zone.  This  is  algebraically 
equivalent  to  solving  the  Rankine-  Hugoniot  shock  relations,  except  that  A 
appears  as  a  parameter  and  it  takes  on  the  values  ,  0  <  A  <  1  as  it  describes 
partial  reaction. 

Once  again  one  can  solve  this  problem  graphically  in  the  p,  u-plane,  (u  = 
1/p)  i.e.  by  finding  the  intersection  of  the  Rayleigh  line  and  the  partial- 
reacted  Hugoniot.  As  before  the  Rayleigh  line  is  given  by 


and  the  Hugoniot  is  given  by 

H:  e{p,v,X)  +  pv+^{v/vQ)^D^  =  (65) 

which  for  the  ideal  EOS  becomes 
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The  negative  of  the  slope  of  the  Rayleigh  line  is  proportional  to  the  square 
of  the  detonation  speed,  D"^. 

If  one  plots  the  Hugoniot  in  the  p,u  plane  for  different  values  of  A,  the 
following  points  hold  for  most  standard  explosive  equations  of  state. 

•  For  partial  to  complete  exothermic  reaction,  0  <  A  <  1,  the  Hugoniot  is 
displaced  upwards,  away  from  the  origin,  from  the  shock  Hugoniot(H; 
evaluated  with  A  =  0.) 

•  The  intersection  of  a  Rayleigh  line  and  the  complete  reaction  Hugoniot, 
(H:  evaluated  with  A  =  1),  generally  defines  two  solutions  on  the  shock 
branch.  The  high  pressure  solution  is  called  the  Strong  end  state,  the 
lower  pressure  solution  is  called  the  Weak  end  state. 

•  At  a  sufficiently  small  value  of  the  detonation  speed,  or  the  Rayleigh 
line  and  the  Hugoniot  are  tangent  and  the  two  values  for  the  end  states 
merge.  There  are  no  solutions  for  values  of  D  below  this  critical  value. 
The  minimum  detonation  speed  is  called  the  Chapman-Jouguet  deto¬ 
nation  velocity,  with  D  =  Dcj. 

Figure  6.  shows  a  sketch  of  the  partial-reaction  Hugoniots,  Weak  and 
Strong  and  CJ  end  states  in  the  p,u-plane  for  different  detonation  speeds. 
Note  that  since  the  Rayleigh  line  is  derived  simply  from  the  conditions  of 
mass  and  momentum  conservation  throughout  the  detonation  structure,  it 
holds  for  each  value  of  A.  The  detonation  states  must  lie  on  the  Rayleigh 
line. 

For  a  given  value  of  D,  such  that  the  Rayleigh  intersects  the  complete  re¬ 
action  Hugoniot,  the  detonation  structure  can  be  described  as  follows.  There 
are  two  possible  structures.  The  Weak  structure  has  no  shock  and  starts 
from  the  initial  state,  point  0  in  figure  6,  and  goes  up  in  pressure  along  the 
Rayleigh  line  to  the  Weak  state.  Although  this  structure  is  admissible  from 
the  point  of  view  of  being  a  valid  solution  to  the  governing  equations,  the 
standard  argument  used  to  rule  it  out  is  that  there  is  no  precursor  ignition 
mechanism  available  to  ignite  the  reaction  zone. 
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Figure  6:  The  p,t;-plane  showing  partial  reaction  Hugoniot  and  the  Weak, 
Strong  and  CJ  state  admitted  by  the  R-H  algebra.  Rayleigh  lines  are  shown 
as  broken. 

The  Strong  structure  is  called  the  ZND  solution,  and  starts  out  with  an 
inert  shock.  Thus  in  the  p,u  plane,  the  state,  jumps  from  point  0,  to  the  inert 
shock  point  at  point  N,  at  the  intersection  of  the  inert  Hugoniot  (A  =  0)  and 
the  Rayleigh  line.  For  realistic  equations  of  state,  the  shock  pressure  is  the 
highest  pressure  within  the  detonation  structure.  The  shock  also  provides 
the  mechanism  for  the  ignition  of  the  reaction.  The  pressure  then  drops, 
from  its  highest  value  at  state  N,  to  its  termination  value  at  the  Strong  state 
at  point  S  as  A  increases  from  0  to  1. 

For  a  simple  equation  of  state,  these  calculations  can  be  done  analytically, 
however  for  a  more  complex  EOS  the  solution  to  the  algebra  must  be  carried 
out  numerically.  Recall  that  the  solution  for  the  end  state  pressure  was  in 
fact  given,  in  an  earlier  section,  for  the  ideal  EOS  by  equation  (25).  Notice 
that  the  Strong  and  Weak  state  coalesce  to  the  CJ  state  when  the  argument 
of  the  square  root  vanishes. 

4.2  Exercise:  Calculation  of  Dcj  for  the  ideal  EOS 

Show  that  the  condition  that  determines  the  CJ  detonation  velocity  is 
easily  obtained  in  the  strong  shock  limit  as 


24 


^CJ  ~  2(7^  ~ 


(67) 


4.3  ZND  Spatial  Structure 

To  find  the  spatial  structure  of  the  detonation,  one  still  has  to  integrate  the 
rate  equation.  We  use  a  *  subscript  to  denote  the  solution  of  the  partial 
reaction  Rankine-Hugoniot  algebra  for  functions  of  A 


p  =  p.(A),f/  =  t/.(A),p  =  p.(A). 


Then  the  rate  equation  is  written  as 


(68) 


^  _  r{p„{X),p,{X),X) 
dn  ~  (/.(A) 


which  in  turn  can  be  integrated  to  obtain  the  spatial  distribution  n  =  n(A) 
as 


u.a) 


(70) 


4.3.1  Exercise:  Calculation  of  the  ZND  reaction  zone  structure 
for  an  ideal  EOS 


Show  that  for  the  ideal  equation  of  state  that  the  CJ,  ZND  solution  structure 
can  be  described  by  the  relations 


where  i  is  defined  by 


7  +  1 


P.  =  Po- 
f/.  =  -Dcj 


7-^ 

7  +  1’ 


p,  —  PqDqj 


ill 

7  +  1’ 


(71) 

(72) 

(73) 


(74) 
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and  that  the  distribution  of  the  reactant  is  found  from  the  integral  of  the 
rate  equation 


Carry  out  the  integration  for  the  case  when  the  rate  law  is  given  by  sim¬ 
ple  depletion  alone  and  the  premultiplier  of  the  rate  constant  has  no  state 
dependence.  Then  r  is  of  the  form 

r  =  fc(l  —  A)",  with,  k^v  =  constant.  (76) 

5  Intrinsic  Geometry  and  Shock- Attached  Co¬ 
ordinates 

We  next  discuss  the  use  of  intrinsic,  shock-attached  coordinates,  in  order 
to  describe  curved,  time-evolving  detonation  waves.  Ultimately  we  must 
represent  the  wave  system  in  a  fixed  laboratory  frame  which  we  take  as 
{x,y,z).  The  intrinsic  coordinate  frame  we  discuss  here  has  the  following 
properties: 

V 

•  The  curvilinear  coordinates  are  normal  to  the  shock  and  parallel  to  the 
shock  surface. 

•  One  of  the  coordinate  surfaces  is  coincident  with  the  shock  surface. 

•  The  reaction  zone  follows  immediately  behind  the  shock  and  is  ac¬ 
counted  for  in  these  calculations. 

•  The  coordinates  reflect  the  growth  (or  shrinkage)  of  the  shock  area  (or 
the  arc-length)  in  an  expanding  (converging)  geometry. 

•  The  coordinates  are  more  complicated  in  their  complete  form,  than 
Cartesian  coordinates  (say).  This  can  be  a  disadvantage. 

There  are  two  types  of  intrinsic  coordinates  that  have  been  developed 
extensively  for  the  description  of  Detonation  Shock  Dynamics;  a  2D  version 


26 


Figure  7:  A  sketch  of  arc-length  angle  coordinates  used  in  2D  for  DSD 

that  uses  arc  length  along  the  shock  and  the  angle  of  the  shock  normal  rela¬ 
tive  to  a  fixed  direction,  and  a  3D  version  that  uses  the  orthogonal  coordinate 
net  defined  by  surfaces  parallel  to  the  shock  surface  and  the  lines  of  prin¬ 
ciple  .curvature  defined  in  the  shock  surface  itself.  The  first  coordinates  are 
described  completely  in  a  recent  report  by  Bdzil  and  Fickett  [9].  The  sec¬ 
ond  will  be  described  in  these  notes.  The  presentation  will  be  specialized,  in 
places  to  2D,  to  shorten  the  presentation.  The  additions  required  for  3D  are 
straightforward.  A  representative  sketch  of  the  coordinate  systems  is  shown 
in  figures  7  and  8  representative  of  the  2D  and  3D  coordinate  systems  used 
in  the  theory. 

5.1  Bertrand  coordinates 

The  coordinates  that  we  pick  are  known  as  Bertrand  coordinates.  The  start¬ 
ing  point  is  a  shock  surface  that  can  be  represented  quite  generally  in  terms 
of  laboratory-fixed  coordinates  (i,  y,  z)  as 

t/»(x,j/,z;t)  =  0.  (77) 

This  equation  constrains  the  lab-coordinate  position  vectors  in  the  surface 
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Figure  8:  A  sketch  of  the  orthogonal  net  comprised  of  surfaces  parallel  to 
the  shock  and  the  surfaces  aligned  with  the  directions  of  principle  normal 
curvatures  in  the  shock  surface,  used  for  3D,  DSD  theory 

to 


x  =  Xs.  (78) 

The  normal  to  the  surface  can  always  be  calculated  by  the  formula  (the  sign 
is  chosen  to  be  in  the  direction  of  the  unreacted  explosive) 


ivv»r 


(79) 


The  shock  surface  can  be  represented  in  terms  of  any  valid  surface  pa¬ 
rameterization  that  we  choose,  and  we  suppose  that  the  shock  surface  can 
be  represented  by  f  =  Xa{^\,^2yt)y  The  intrinsic  coordinates  are  related  to 
the  laboratory  coordinates  by  the  change  of  variable  given  by 


X  =  Xs(6,6,<)  +  ««(6,6;<)>  (®^) 

where  the  variables  n,^i,42  are  respectively  the  distance  measured  in  the 
direction  of  the  normal  to  the  shock  wave,  and  the  arclength  measured  in  the 
directions  defined  by  the  principle  normal  curvatures  at  the  shock  surface. 
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If  we  travel  in  the  shock  surface  in  the  direction  of  or  ^2)  we  define  the 
unit  vectors  of  the  local  basis  of  this  coordinate  system  in  the  shock  surface 
(which  are  the  tangent  vectors  to  the  coordinate  directions),  namely 

1  —  X  _ 

We  restrict  the  remaining  development  to  2D  to  simplify  the  presentation. 


5.2  Arc  -  length  angle  coordinates 

Let  (t>  measure  an  angle  of  the  shock-attached  normal  to  a  fixed  direction, 
here  aligned  in  the  z-direction  (say).  In  the  most  general  case,  the  reference 
direction  can  be  a  space  curve  with  arbitrary  orientation.  We  will  refer 
to  the  reference  space  curve  as  the  edge.  The  arclength  in  the  surface  is 
measured  relative  to  the  intersection  of  the  shock  surface  with  the  reference 
curve,  hence  this  intersection  is  the  instantaneous  origin  for  the  intrinsic 
coordinate  system.  This  definition  of  the  origin  allows  for  a  unique  mapping 
of  the  coordinates  of  the  intrinsic  coordinates  to  laboratory  coordinates.  In 
particular  n  =  0  always  defines  the  shock  surface,  and  ^  =  constant  are  space- 
time  trajectories  that  measure  constant  arclength  on  the  shock,  relative  to 
the  edge. 

At  each  point  on  the  shock,  the  angle  ^  is  defined  so  that  the  function 
<j>{i,t)  describes  the  orientation  of  the  normals  of  the  2D  shock.  The  normal 
and  tangent  unit  vectors  are  defined  in  terms  of  the  lab-coordinate  basis 
vectors  Cj,,  in  the  x  and  z-directions  by 


n  =  sin{<f>)ex  +  cos{<l>)ez,  i  =  cos(^)ej;  —  sm(^)e2.  (82) 

The  Frenet  formulas  applied  to  this  coordinate  system  relate  the  derivatives 
of  the  unit  vectors  along  the  coordinate  directions  and  reflect  the  fact  that 
our  coordinate  system  is  locally  orthogonal  and  is  aligned  along  the  lines  of 
principle  normal  curvature.  These  formulas  in  2D  are 


.  dn  f 
-«n,  —  =  Kt. 
OK 


(83) 


Notice  that  in  2D,  the  derivative  of  the  angle  <f>,  with  respect  ^o  the  arclength 
^  defines  the  curvature  k 
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(84) 


d<t> 


As  an  exercise  the  reader  can  check  that  according  to  the  above  formulas  the 
basis  of  this  intrinsic  coordinate  system  is  locally  orthogonal,  (i.e.  the  dot 
product  of  the  normal  and  tangent  vectors  is  zero). 

The  relation  between  the  laboratory  coordinate  x  and  the  arclength,  angle 
coordinates  is  as  follows.  Since  positions  on  the  shock  are  described  by 

dx  =  id(,  (85) 

then  it  follows  from  the  formula  for  i  that 

dx  —  cos{(j))d(,  dz  =  — sm(^)d^,  (86) 

and  integration  of  these  formulas  then  obtains 

X  =  Xe(t)  +  cos(<^(i,  i))d^,  z  =  Ze(t)  -  t))d^.  (87) 

The  coordinates  of  the  edge,  Xe{t),Ze{t),  are  introduced  for  the  purpose  of 
measuring  the  arclength  from  a  unique  origin  in  the  surface  at  each  instant 
of  time.  These  coordinates  are  found  by  locating  the  intersection  of  the  edge 
and. the  shock  surface  and  thus  are  necessarily  functions  of  time.  In  general, 
the  edge  coordinates  must  be  determined  in  the  course  of  solving  for  the 
evolution  of  the  shock. 

5.3  Change  of  variable 

The  coordinate  transformation  implied  by  the  definitions  of  the  previous 
section  change  the  representation  of  the  governing  equations  that  we  have  to 
solve.  A  change  of  variable  is  required  to  transform  the  equations  into  the 
n,^,f  coordinates.  In  what  follows,  when  a  derivative  is  for  fixed  laboratory 
position,  we  will  used  the  notation  jj,  and  for  fixed  position  in  the  shock 
attached  coordinates,  we  will  use  In  particular  we  must  show  how 

Cii  Cz  goes  to  i,  n,  (88) 
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d  d  '  d  d 

dx'  dz  ^  dn'  di  ' 

and 

(89) 

d  ^  d 

For  example  we  need  to  compute  the  V  operator,  which  in 
coordinates  is  written  as 

(90) 

the  lab-fixed 

By  the  chain  rule  we  have 

(91) 

d  di  d  dn  d 

dx  dx  d(  "**  dxdn 

d  d^  d  dn  d 

dz  dzdi^  dzdn 

(92) 

and  we  compute  the  derivatives  from  the  coordinate  transfor¬ 

mation. 

tf  we  differentiate  the  coordinate  transformation  (in  2D) 

xex  + ze^  =  x,{(,t)  +  nh{^,t),  (93) 

with  respect  to  x  and  z  respectively,  and  use  the  definition  of  i  and  the  Frenet 
formulas  (83),  we  obtain  the  following  formulas 

id(dn  d(f 

(94) 

;9(,9n  d(. 

'■  =  ‘a;  +  il"  + 

(95) 

If  we  take  the  dot  product  of  the  above  formulas  with  n  and  i  respectively, 
we  get  the  following  relations  for  the  direction  cosines  between  the  lab-fixed 
coordinates  and  the  intrinsic  coordinates 
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*  .  dn  . 

e,.n  =  ^,  e, 

.  *  dn  . 

e.-n  =  -,  e. 


.(  =  (l+nK)|i, 
■  f  =  (1  +  nK)|i. 


(96) 


The  lab-fixed  unit  vectors,  e^,  Cj  can  be  written  in  terms  of  their  components 
in  the  intrinsic  frame  as 


e®  =  (CiE  '  n)n  -f  (cj,  •  i)i, 

Cz  =  (cj  •  n)n  -f  {tz  •  i)i.  (97) 

The  V  operator  can  then  be  written  as 


V  = 


[(e..n)n  +  (e..J)(](gA  +  ^|.j 
[(e,  •  n)n  -h  (cz  •  i)i] 


+ 


dzd(  dz  dnj 


(98) 


If  we  multiply  out  the  above  result  and  use  the  properties  of  the  unit  vectors 
thaf  |f|  =  ln|  =  I,  <  •  n  =  0,  or 


and  use  the  fact  that  dijdn  =  0,  then  the  following  formula  results 

^  + "r-- 

1  -|-  UK  on 

As  another  example,  by  using  the  definition  of  the  velocity  in  the  intrinsic 
coordinates 


u  =  U(i  -I-  u„n. 


(101) 
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we  calculate  the  divergence  V  •  u,  as 


V  •  u  = 


1  du^  ^  «n 
1  +  n/c  54  5n  1  +  UK 


—  «Uf. 


(102) 


The  time  derivative  in  the  lab-fixed  coordinates  is  related  to  that  in  the 
shock-attached  coordinates  by  the  formula 


£|  _£|  4.^1  A  +  A 


(103) 


where  dnldt\g  =  —Dn  can  be  identified  as  the  negative  of  the  normal  com¬ 
ponent  of  the  shock  surface  velocity  and,  d(/dt\g  =  B  is  the  instantaneous 
rate  of  increase  or  decrease  of  arclength  along  the  shock.  Thus  we  write 


_  A|  _  n  A+fiA 


(104) 


5.4  The  governing  equations  in  intrinsic  coordinates 

The  governing  equations  in  intrinsic  coordinates  follow  from  a  straightforward 
application  of  the  formulas  of  the  previous  section.  We  summarize  them 
below:  The  equation  of  conservation  of  mass  is  given  by 


1 


+  nK  1  -f  nK  cf^ 


The  n-momentum  equations  is  given  by 


K 


+  (un  ~  Ai)“x — I-  T"; - ^CTT 

dt  an  I  +nK  d(  ^  1  -h  uk 


-  Bu(K 


dp  - 


The  ^-momentum  equation  is  given  by 


(106) 


du 


du^  U(  du^ 


dt  an  1  +  UK  di 


+  lt«Wn 


1  4-  UK 

1  dp 


\  +  UK  d^ 


=  0. 


(107) 
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The  energy  equation  is  given  by 


dt  ,  r.  uede  de 

+  (Un  —  ^n)-^  +  7"7 — 

at  an  1  +  nK  a?  a^ 


I  rn  n  I  1 

"oT  "h  v“n  ~  ^nj'5 - T  T~;  ^ 

at  an  1  +  n/c  a^ 


=  0. 


The  rate  equation  is  given  by 
dX 


dX 


u(  dX 


.dX 


^  +  (u„-D„)^+-^^  +  S^  =  r. 
at  an  I  +  nn  a^ 


(108) 


(109) 


5.5  Additional  kinematics  of  the  surface 

In  the  equations  of  the  last  section  the  quantity 

B-|W,  (110) 

appears  and  represents  the  instantaneous  change  in  the  arclength  along  the 
shock.  Here  we  give  the  derivation  to  show  how  this  term  is  described  com¬ 
pletely  in  terms  of  quantities  defined  in  the  surface.  By  doing  so,  we  derive 
thej:tnema<tc  surface  relations. 

Differentiating  the  change  of  variable  formula 


xt  =  Xsitt) +nh{(,t). 


(111) 


with  respect  to  t,  holding  xt  fixed  (and  using  the  chain  rule  and  Frenet 
formulas)  obtains, 


0  = 


dxs.  ,dn 


(112) 


We  specialize  this  result  to  the  surface  x  =  x,,  set  n  =  0,  and  use  the 
definition  dx^jd^  —  i  and  dnfdt\£=  —Dn  to  obtain 


Next  we  differentiate  (113)  with  respect  to  holding  t  fixed,  and  use  the 
Frenet  formulas  to  obtain  the  following  vector  equation  that  describes  the 
kinematics  in  the  shock  surface 


dPn 


h  —  DnKt  =  0. 


(114) 


We  dot  the  last  equation  with  n  to  obtain  the  n-  component  and  obtain 


n 


di. 


^1  _^-n 
at'*- 


(115) 


and  dot  the  same  equation  by  i  to  obtain  the  (  -  component  as 


This  pair  of  surface  relations  can  be  used  to  describe  the  shock  surface 
evolution  solely  in  terms  of  angle  -  arclength  coordinates  as  follows.  From 
the  formulas  that  define  the  normal  and  tangent  vector  in  terms  of  the  shock 
normal  reference  angle  (82),  and  the  definition  of  the  plane  curvature, 
(84),  we  find  that 


di.  M 

#■=-"¥ 

By  integrating  equation  (116)  with  respect  to  and  using  the  definition 
K  =  d<f>ld^,  we  get  the  explicit  expression  for  the  instantaneous  rate  of  change 
of  the  arc  length. 


where  f{t)  is  an  integration  constant.  Finally  from  the  Dn  —  k  relation,  we 
have  the  fact  that  Dn{ii)  is  presumed  to  be  a  known  function.  Substituting 
these  results  back  into  (116)  yields  the  intrinsic,  integro-differential  equation, 
for  <!>((,. t). 

di  dPM  dig. 

at  d(  die  dP'  '  ' 

This  equation  was  derived  first  in  [16]  by  a  different  route. 
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Importantly,  the  results  of  this  section  can  be  used  to  make  explicit  es¬ 
timates  of  the  asymptotic  orders  of  certain  dynanoic  terms  that  occur  in  the 
governing  equations.  We  explore  this  in  the  next  section. 

6  Asymptotic  Scaling 

Recall  the  basic  results  that  we  had  previously  obtained  by  discussing  Carte¬ 
sian  coordinates  where  the  shock  locus  was  given  by 


=  z  —  Dcjt  -  Zs(x,  t)  =  0.  (120) 

The  normal  detonation  velocity  was  given  by 


_  Dcj  +dzgldt 

[I +  {dzjdxyY/^^ 

and  the  plane  curvature  was  given  by 


(121) 


d^Za  1 

dx^  [1  -f  (dz,/ax)2]3/2- 


(122) 


Also  recall  that  the  steady  ZND,  1/2-reaction  zone  length  was  given  by  the 
integral 


U.jX)  , 

A=o  r(p.(A,p.(A),A)  ’ 

where  the  *  subscript  refers  to  quantities  defined  by  the  steady  state. 

In  what  follows,  within  this  section,  we  adopt  the  notation  convention 
where  a  quantity  with  a  ()  refers  to  a  dimensional  quantity  and  the  quantities 
without  a  tilde  are  dimensionless  quantities  that  are  scaled  with  respect  to 
the  dimensional  unit.  In  particular  we  choose  the  following  scales  for  length, 
velocity  and  time 


—  Dcj, 

—  hzIDcj- 


length 

velocity 

time 
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(124) 


Then  the  formulas  (in  Cartesian  coordinates)  for  the  normal  detonation  ve¬ 
locity  and  the  curvature  become 


Pn  \+dz,ldt 

bcj  [i  +  idzjaxni^' 

and 


d^zjdx'^ 

[1  -f  (dz,/aa:)2]3/2- 


Recall  the  limit  of  small  shpck  slope,  can  be  expressed  as 


(125) 


(126) 


^«0(e).  (127) 

The  reaction  zone  thickness  is  now  the  basic,  0(1),  unit  of  distance  measure¬ 
ment  in  our  theory.  The  natural  way  to  express  slow  variations  along  the 
shock,  measured  in  a  direction  transverse  to  the  shock,  is  to  introduce  the 
slowly  varying  length  scale 


X  =  ex.  (128) 

Transverse  variations  are  measured  on  the  X  0(1)  length  scale  so  that 
when  djdX  ~  0(1)  then  djdx  —  edfdX  ~  0(c). 

It  follows  that  the  shock  slope  scales  as 


dz,  ^  dzt 
dx  dX 


(129) 


The  argument  that  we  had  used  in  section  3.5,  suggests  that  if  temporal 
changes  in  the  movement  of  the  detonation  shock  are  to  affect  the  detonation 
velocity  at  the  same  order  as  the  shape  changes,  then  the  relevant  time  scale 
must  be  such  that 


dz, 

dt 


( dz 

O(c^),  since  ( j  ~  O(e^).  (130) 


Thus  we  introduce  the  slow  time  variations  and  introduce  the  time  scale 


r  =  c^t. 


(131) 
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so  that 


2± 

di~‘  ar 


(132) 


6.1  The  expansion  procedure 

What  follows  next  is  a  discussion  of  the  expansion  procedure  that  is  used  to 
derive  a  reduced  set  of  equations  that  will  describe  the  quasi-steady,  weakly 
curved  dynamics  of  near-CJ  detonation  shocks.  The  smallness  of  the  shock 
curvature  is  measured  relative  to  the  to  the  reaction  zone  thickness  in  the 
normal  direction. 

There  are  two  main  developments  that  must  be  in  place  before  we  can 
start  this  reduction,  and  they  have  been  covered  in  the  previous  sections:  i) 
The  governing  equations  need  to  be  written  in  a  shock-attached  frame,  ii) 
Scaling  assumptions  that  describe  the  nature  of  the  unsteadiness  and  their 
relation  to  the  shape  changes  of  the  shock  must  be  adopted.  The  scaling  as¬ 
sumptions  are  those  associated  with  a  particular  asymptotic  distinguished 
limit  of  the  reaetive,  multidimensional  Euler  equations.  The  multi-scale 
asymptotic  expansions  that  are  described  next  (although  perhaps  robust  in 
their  ability  to  model  the  underlying  physical  system)  are  mathematically 
restricted  to  describe  the  shock  dynamics  within  the  limits  of  the  scaling 
assumptions  and  following  expansions. 

The  expansions  that  we  list  here  are  motivated  by  the  dynamics  implicit 
in  the  shock  relations  and  the  fact  that  to  leading  order,  we  expect  the 
detonation  to  behave  in  a  quasi-steady  fashion.  We  reintroduce  dimensional 
equations  here,  and  temporarily  retain  c  as  the  small  parameter  that  describes 
the  small  relative  curvature.  The  variable  n  still  measures  distance  in  the 
shock  normal  direction  and  is  negative  within  the  reaction  zone  and  zero  at 
the  shock.  We  retain  t  as  the  relevant  time  scale,  and  use  C  slow 

transverse  spatial  scale  that  measures  arclength  along  the  shock,  defined  by 


C  =  (133) 

The  expansions  of  the  dependent  variables  for  the  simplest  version  of 
DSD-theory  are 


=  eu[^^(n)  -f  o(c). 
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(134) 


p  =  /)W(n)  +  cV^^(n,C,r)  +  o(c^), 

p  -  pW(n)  +  eV*H«,C»T)  +  o(€^)> 

A  =  A^®^(n,  t)  +  c^A^®^(n,  (,  r)  +  o{c^) 

In  addition  to  the  dependent  variable,  we  also  expand  the  normal  detonation 
velocity  (in  the  simplest  version  of  the  theory)  in  terms  of  the  shock  curvature 
as 


£)„  =  l+e^DW(C,r)  +  ...,  (135) 

where  k  is  related  to  e  by 

K  =  (136) 

Note  by  assuming  the  expansions  (135)  and  (136)  we  have  assumed  that 
the  leading  order  correction  to  the  normal  detonation  velocity  is  0{k).  This 
assumption  is  verified  if  the  analysis  determines  the  coefficient  D^\  and  is 
correct  for  the  simplest  case  presented  here.  However  in  general  the  leading 
order  correction  to  Dn  can  have  a  different  dependence  on  k.  For  example, 
in  [5]  and  [8]  it  is  shown  that  the  leading  correction  to  Z)„  is  0{k  In  «)  for  a 
rate  of  the  form  r  =  k{p,p){\  —  A). 

The  0(1)  terms  in  the  expansions  for  Un,p,p  and  A,  reflect  the  steady, 
ID,  ZND,  solution  in  the  shock-attached  coordinates.  The  full  expansions 
for  the  dependent  variables,  the  curvature  and  the  detonation  velocity  are 
substituted  into  the  governing  equations  (105)  -  (109).  This  leads  to,  in 
the  standard  way,  a  hierarchy  of  systems  of  equations  at  0(1),  0(c),  and 
O(e^).  E2ich  set  of  equations  must  then  be  solved  in  a  consistent  manner 
to  ultimately  produce  a  uniform  approximation  within  the  space  and  time 
domain  of  interest. 

While  this  procedure  is  general  and  precise,  we  replace  it  with  a  sim¬ 
ple  and  entirely  equivalent  procedure,  whereby  we  analyze  each  term  in  the 
governing  equations,  using  estimates  in  terms  of  c  that  are  found  in  the  ex¬ 
pansions  (136)  and  from  the  scaling  definitions.  We  then  discard  terms  that 
are  smaller  than  O(e^),  i.e.  below  the  order  of  the  curvature.  This  analysis 
was  carried  out  systematically  in  [8],  however  Whitham’s  shock  ray  coordi¬ 
nates  are  used  rather  than  the  Bertrand  coordinates  that  appear  in  these 
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notes.  What  follows  is  a  discussion  that  analysis  and  we  go  through  the  gov¬ 
erning  equations,  essentially  term  by  term.  We  discuss  the  negligible  terms 
in  detail. 

Starting  with  the  unsteady  term  in  (105),  we  get  the  estimate 

(137) 

This  estimate  follows  from  two  facts,  that  the  leading  order  approximation 
is  assumed  not  to  depend  on  r  and  that  the  next  correction  to  p  is  O(e^). 
The  term 


(138) 

1  -f  n/c 

because  we  have  slow  transverse  variations  and  is  assumed  to  be  0(c)  and 
not  to  depend  explicitly  on  C.  The  term 


KUf  =  O(c^).  (139) 

The  term  in  the  mziss  conservation  equation  is 

=  (140) 

and  to  show  this,  we  need  to  invoke  an  estimate  of  B  =  d^/dty. 

Recall  from  (116)  that  B  must  satisfy  the  kinematic  surface  relation 

^  =  (141) 

OK 

so  that  in  order  to  be  consistent, 

S=^|,~0(e).  (142) 

A  list  of  other  negligible  terms  and  their  estimates  obtained  by  similar 
arguments  is  found  below.  In  the  n-momentum  equation  we  find 


u|«: 

1  -|-  n/c 


dun. 

dt 


0(6“), 


dUn 
1  -b  n/c  dK, 


0(e“), 


0(€'‘),  Bu(K  = 


(143) 
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The  estimates  of  negligible  terms  in  the  momentum  equation  are 


dt 


U4  du( 
1  +  UK  d( 


0(£*), 


=  O(c’),  =  0(e»),  Bu,«  =  0{^). 

1  +  n/c 


In  the  rate  equation  we  neglect 


(144) 


From  the  ^-momentum  equation  we  find  the  following,  simple  result. 
From  the  0(c)  equation  we  find  that  obeys 


¥  =  «■ 


(146) 


Integrating  this  equation  shows, 


^"  =  /({,r).  (147) 

But  since  the  normal  shock  relations  show  that  there  is  no  jump  in  across 
the ‘shock,  combined  with  the  assumption  that  the  upstream  state  is  quies¬ 
cent,  shows  that 


=  0,  for  all  (  and  t.  (148) 

Consequently  we  find  that 


U5  =  o(c),  (149) 

and  some  of  the  estimates  of  this  section  that  involve  can  be  further* 
refined. 


7  The  Reduced  Problem  for  Dn  —  k, 

If  we  retain  only  the  terms  in  (105)  -  (109)  that  may  contribute  terms  up  to 
order  O(c^)  ~  we  find  the  reduced  set 
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Q 

— (/)(Un  -  D„))  +  pUnK  =  0, 

(ISO) 

.  _  .dun  dp  ^ 

/)(«n-Dn)— +  —  -0, 

(151) 

dn  p^  dn  ' 

(152) 

/  r,  ^9X 

(153) 

In  addition  to  these  equation  we  add  the  (strong)  shock  relations  at  the 
location  n  =  0,  which  we  reproduce  here  for  convenience 


P>  =  P 

Un  =  Un-  Dn  = 


7  +  1 

7-1’ 
7-1 
7  +  1’ 


Ps  = 


-PoDl, 


7  +  1 
=  0,  A,  =  0. 


(154) 


The  shock  relations  essentially  serve  as  boundary  conditions  for  fixed  Dn. 

As  part  of  the  general  formulation,  a  total  integral  of  the  energy  can  be 
obtained  from  combinations  of  the  above  differential  equations.  This  relation 
is  essentially  Bernoulli’s  equation  for  the  total  energy  on  a  streamtube,  where 
the  value  of  the  total  energy  integral  is  evaluated  at  the  shock.  The  strong 
shock  result  is  given  by 


^{PiPi^)  +  “  + 

P 


(155) 


7.1  Derivation  of  the  master  equation 

In  this  section,  we  give  a  derivation  of  the  master  equation  for  the  general 
e(p,  u.  A)  equation  of  state,  and  then  take  the  special  case  of  the  polytropic 
equation  of  state.  We  start  with  an  arbitrary  EOS  to  obtain 


de  dp  de  dp  de  dX 
dp  dn  dp  dn  dX  dn 


i.^  =  0 

dn  ’ 


(156) 
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Which  we  rewrite  as 


1p/p^  -  aejdA  dp  dp  dt/dx  ax 

dejdp  dn  dn  defdpdn 

We  note  the  general  definition  of  the  frozen  sound  speed  as 

2  ^  IpIp"^ -del dp] 

""  -  dt/dp  ’ 

so  that  the  previous  result  can  also  be  written  as 


(157) 


(158) 


^  de/dX  / _ 
dn  dn  dejdp  \^/ 


(159) 


The  next  step  is  to  eliminate  the  pressure  p  and  the  density  p  in  favor 
of  the  relative  velocity  (/„  =  «n  -  Dn  and  the  mass  fraction  A.  From  the 
momentum  equation  we  find  that 


dp  ..  dUn 
dH  "  dn  ’ 

and  from  the  mass  conservation  equation  we  find  that 


(160) 


dp  Un  +  Dn  P  dUn 

dn  ~  U,  l/„  dn  ■ 


(161) 


These  two  substitutions  are  then  used  in  the  energy  equation  and  results  in 
an  equation  that  is  quasi-linear  in  Un  and  A.  If  we  use  A  as  independent 
variable  instead  of  n,  the  rate  equation  provides  the  replacement 


d  _  r  d 


(162) 


By  replacing  the  derivative  djdn  in  the  quasi-linear  form  of  the  energy 
equation  by  djdX  and  replacing  dXIdn  by  r/Un  and  then  solving  for  dUn/ dX, 
we  obtain  an  equation  that  determines  t/„  as  a  function  of  A 


dUn  _  p/>(^n  +  £>n)«+7||^^ 

dX  ~  r[c2-(/2] 


(163) 
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7.2  The  special  case  of  the  ideal  EOS 

Now  we  turn  to  an  example  of  the  ideal  equation  of  state.  This  choice  has 
the  advantage  that  much  of  the  theory  has  a  simpler  form,  and  that  there  is 
a  large  body  of  established  theory  and  niimerics.  The  ideal  (or  polytropic) 

EOS  is  given  by 

P  1 

e  =  ^ - --QX. 

P7-  1 

(164) 

Then  some  of  the  thermodynamic  quantities  in  the  expressions  of  the  previous 
section  are 

Ide/aX 

pde/dp 

(165) 

P 

(166) 

By  using  the  total  energy  equation  we  can  get  an  expression  for 
depends  on  (/„  and  A  as 

that  only 

(167) 

The  master  equation  can  then  be  written  in  the  form 

at/n  _ 

dX  rr}' 

(168) 

where 

=  (7-l)gr-c=*(f/„  +  £>n)/c, 

(169) 

and 

(170) 

* 

A  rate  law  with  simple  depletion  that  has  a  premultiplying  rate  constant  that 
is  sound  speed  (temperature)  dependent  is  generally  of  the  form 

r  =  ik(c2)(l  -  A)*'. 

(171) 

.. 
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7.3  The  shooting  problem  for  the  Dn  —  k,  relation. 

Here  we  discuss  the  shooting  problem  that  determines  the  Dn  —  n  relation. 
Note  that  the  sign  convention  we  have  used  assumes  that  when  the  detonation 
shock  is  convex,  relative  to  the  ambient  flow,  that  the  curvature  is  positive. 
This  is  the  case  of  a  diverging  detonation. 

The  master  equation  (163)  was  derived  for  an  arbitrary  EOS  of  the  form 
e(p,  p.  A),  and  has  the  general  form 


where 


dU.  _  t/n^ 

dX  rrj  ' 


(172) 


♦  =  -i^^-C=((/n  +  On)<=,  (173) 

p atlap 

and  Tj  is  defined  by  (170).  The  solution  to  the  master  equation  can  be  entirely 
described  in  a  (/„  —  A  plane  if  the  right  hand  side  of  (172)  can  be  written 
entirely  as  a  function  of  at  most  Un  and  A.  If  this  is  not  the  case,  and  for 
example,  the  density  dependence  cannot  be  eliminated  explicitly  from  (172), 
then  the  continuity  equation  must  be  added  for  a  complete  description. 

For  the  shooting  problem,  for  general  EOS,  one  must  also  solve  the  shock 
rela,tions  to  obtain  a  boundary  condition  of  the  form 


Un  =  Un{Dn)  at  A  =  0.  (174) 

The  integral  curve  must  start  at  the  shock,  with  A  =  0  and  for  positive 
curvature  pass  through  a  singular  point  that  corresponds  to  a  condition  where 
the  flow  becomes  sonic  (i.e.  1/  =  0).  So  that  we  have  a  nonsingular  solution  at 
the  sonic  point,  the  integral  curve  must  be  such  that  $  and  77  simultaneously 
vanishes  there.  The  condition  that  the  Un  -  X  integral  curve  pass  through 
the  point  in  the  Un  —  X  plane  defined  by  the  conditions 


$  =  0  and  77  =  0,  (175) 

is  called  the  generalized  CJ  condition  after  Wood  and  Kirkwood,  [3],  and 
must  be  applied  in  addition  to  the  shock  boundary  condition.  The  task  of 
solving  for  an  integral  curve  in  the  Un-X  plane  is  overdetermined  unless  one 
allows  for  a  relation  between  Dn  and  /c. 
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The  master  equation  (163)  subject  to  the  shock  boundary  condition  (174) 
and  the  generalized  CJ  condition,  is  a  self-contained  problem,  providing  that 
it  is  possible  to  write  the  thermodynamic  quantities,  c^,  (de/dA)/(/j5e/dp) 
and  the  rate  law  r,  solely  in  terms  of  (/„  and  A,  This  cannot  be  always  be 
done  for  a  non-ideal  equation  of  state  or  rate  law.  Rather  it  is  true  only 
under  some  special  circumstances,  like  the  example  shown  in  section  7.2  for 
the  ideal  equation  of  state,  and  a  rate  law  which  had  a  rate  premultiplier 
that  is  only  function  of  <?. 

Once  the  integral  curve  t7n(A)  is  found,  one  can  recover  the  density  and 
the  pressure  from  the  integration  of 

=  -AV.  +  (176) 

Un  OA 

and 

subject  to  the  appropriate  shock  conditions.  The  spatial  structure  of  the 
reaction  zone  is  then  finally  determined  by 

U 

n  =  /  ^dA.  -  (178) 

^  ■  Jo  r 

7.4  The  shooting  problem  for  general  EOS 

The  general  formulation  of  the  shooting  problem  requires  the  the  simulta¬ 
neous  integration  of  either  the  mass  or  momentum  equations,  in  addition  to 
the  master  equation  to  describe  the  solution  through  the  reaction  zone.  The 
general  formulation  (for  one  reaction  variable)  can  be  written  in  terms  of  two 
coupled  nonautonomous  equations  for  either  {p,P),  {p,[In)oT  {Un,P)-  The 
missing  dependent  variable  is  supplied  through  the  total  energy  integral, 

e(p,P»  A)  +  -  +  -^  =  e(p„pa,0)  + 

p  2  ^  Ps 

The  relevant  equations  for  p,  Un  and  p,  with  A  as  the  independent  variable 
are 
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II 

(Vn  +  Dn)K  +  -]  , 

(180) 

dl/„  _  Un^ 
d\  rr)  ’ 

(181) 

Q:|  03 

II 

(182) 

A  numerical  procedure  is  a  robust  way  to  determine  the  determine  the 
•solution  and  the  Dn  —  k  relation.  A  successful  iterative  procedure  starts 
the  integration  of  the  ODEs  at  the  shock,  and  integrates  toward  the  end  of 
the  reaction  zone.  For  a  fixed  £)„,  ac  pair,  one  of  the  two  of  the  generalized 
CJ  conditions  will  be  satisfied  first,  while  the  other  is  not.  One  uses  the 
residual,  i.e.  the  nonzero  value  of  the  condition,  to  develop  an  iteration  • 
procedure  (like  a  Newton-Raphson  or  secant  method)  to  change  (say) 
systematically.  The  integration  is  repeated  until  the  shock  conditions  and 
the  generalized  CJ  conditions  are  satisfied  simultaneously.  Lee,  Persson  and 
Bdzil  used  this  procedure  recently  to  determine  both  the  —  k  relation  and 
the  rate  law  for  nonideal,  ammonium  nitrate-based  explosive  used  in  mining 
applications,  [13]. 

8  The  Singular  Perturbation  Solution  to  the 
Dn  —  Relation  for  Ideal  EOS 

More  extensive  discussions  of  the  analysis  presented  here  are  found  in  [5], 
and  in  particular  [8]. 

We  reintroduce  the  dimensional  scales  so  that  the  resulting  problem  that 
we  finally  analyze  is  dimensionless.  The  following  quantities  and  parameters 
are  defined  for  the  ideal  equation  of  state  described  in  section  7.2, 

Un  =  Un!  DcJi  c  =  cl  Dcj,  T  =  T^z/ Dcj 

Dn  =  DnIDcj,  K  =  kirz,  ?  =  QIDqj.  (183) 

Having  done  this  we  obtain  the  dimensionless  version  of  the  master  equa¬ 
tion  for  an  ideal  equation  of  state  and  idealized  rate  law  as 
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(184) 


dUn  _  jj  [(7  -  +  ^n)«l 

dX  “  "  r(c2-i/2) 

where  in  terms  of  Un  and  A  is  given  by 

C*  =  ^^(Dl  -  Ul)  +  ,A(7  -  1),  (185) 

with  the  shock  condition 

Un  =  -^^n.  (186) 

7+  1 

To  illustrate  the  theory  we  will  treat  the  simplest  case  where  A:  in  (171)  is  a 
constant  and  is  not  state-dependent. 

The  objective  here  is  to  describe  the  asymptotic  solution  to  this  problem 
in  the  limit  as  k  0.  The  solution  is  an  integral  curve  in  the  (/„  —  A  -  plane 
that  starts  from  the  shock  and  ends  up  at  the  generalized  CJ  point. 


8.1  Location  of  the  critical  point  in  the  C7„  —  A  -  plane 

The  location  of  the  generalized  CJ-point  is  given  by  the  conditions 


$  =  0,  and  T]  =  0. 

(187) 

These  conditions  respectively  are 

{'r-l)qk{l-XY==C^{Un  +  D)K, 

(188) 

and 

II 

(189) 

In  particular,  in  equation  (188),  as  /c  — ♦  0,  we  must  have  A 
we  can  make  the  assumptions  near  the  generalized  CJ  point. 

1,  SO  that 

Dn  =  1  +  +  •••>  A=l  —  A  -|- - 

(190) 

Substitution  of  those  expansion  into  the  previous  expression  gives  approxi¬ 
mate  expressions 
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We  obtain  two  algebraic  conditions  from  this  analysis  near  the  singular  point. 
Notice  that  the  order  of  D'^  is  not  determined  by  this  analysis. 

8.2  Near-shock  expansion 

The  analysis  of  the  structure  of  the  solution  to  this  problem  breaks  into  two 
pieces,  in  the  spirit  of  layer  analysis  of  singular  perturbation  theory.  There 
is  generally  a  near-shock  layer  and  a  layer  near  the  generalized  CJ  point.  In 
these  notes  we  will  treat  the  simplest  case  possible,  when  the  relationship 
between  and  k  is  linear. 

We  suppose  that  the  near-shock  layer  has  a  regular  perturbation  expan¬ 
sion  in  terms  of  k  of  the  form 

(/„  =  -h  . . . ,  (193) 

and  represent  Z)„  by 

Dr,  =  l  +  KD^^^  +  ....  (194) 

The  object  of  the  analysis  is  to  find  a  formula  for  Substitution  of  the 
above  expansions  into  (184)  and  (186)  and  collecting  powers  of  k  gives  the 
problems  for  At  0(1)  we  obtain 

((c‘“l)*  -  -  1)9,  (195) 

subject  to  the  shock  boundary  condition 

t/(0)  =  _IZli  at  A  =  0,  (196) 

7  +  1 

where  (c^°l)^  is  determined  in  terms  of  and  A  from  (185). 
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At  0(k)  we  have 


((,2)U)  _  +  ((c(0)f-(UWf)^  = 

_  (/(o)(c(°))2((/(o)  +  l)/r(°),  (197) 

(where  =  fc(l  —  A)"),  and  is  subject  to 

[/(i)  =  at  A  =  0.  (198) 

Note  that  (c^)^^^  is  also  determined  from  (185)  in  terms  of  and  A, 

and  For  convenience  we  introduce  the  variable 


i  =  VT^.  (199) 

Both  of  the  problems  listed  above  can  be  solved  explicitly. 

The  0(1)  problem  is  exactly  the  same  as  the  ID,  steady,  CJ,  ZND  deto¬ 
nation  wave,  see  section  4.3.1.  The  0(1)  problem  requires  the  solution  of  a 
separable  first  order  ODE  and  the  solution  is  given,  by 


t/(0)  =  _lzi. 

7+1 

The  solution  to  the  0(/c)  equation  for  is 


(200) 


where  F{t)  is  defined  by 


—  f- 

7  +  1  V7 


(201) 


At  this  point  the  leading  order  correction  to  the  D„  —  k  relation  can  be 
derived  by  using  the  principal  of  elimination  of  the  strongest  singularity.  In 
this  case,  a  singularity  exists  in  the  0(k)  expansion  and  must  be  removed  to 
obtain  a  regular  solution  near  the  end  of  the  reaction  zone. 

Note  that  in  the  near-shock  expansion  that  the  solution  for  becomes 
singular  at  the  tail  of  the  reaction  zone  as  A  — 1,  i.e.  as  £  0,  without 
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further  assumption  would  become  asymptotic  to  0(1/^).  Thus  sufficiently 
close  to  the  end  of  the  reaction  zone  the  expansion  becomes  nommiform. 
The  simplest  argument  that  can  be  made  is  that  this  nonuniformity  is  not 
a  characteristic  of  the  physical  solution  and  it  is  not  present  in  the  physical 
flow.  We  can  eliminate  the  nonuniformity  by  choosing  according  to  the 
formula 


£>a)  =  _  r  (203) 

Jo 

With  this  choice  the  solution  for  becomes 

=  -1-  /'  F{i)di  -  .  f  F{i)di 

7  +  Wo  7(7  +  1)  Z/ 

and  this  solution  is  not  singular  as  Z  — >  0  since  ^  F{£)d£  ~  £. 

Notice  also  that  the  integral  that  defines  depends  on  the  form  of  the 
rate  law  and  the  equation  of  state.  For  this  choice  of  the  rate  law  made  here, 
the  integral  is  convergent  only  for  values  of  0  <  1/  <  1.  For  1/  =  1,  the  integral 
is  formally  divergent,  which  indicates  that  the  leading  order  correction  to  Dn 
is  not  0(/c),  but  is  larger.  (It  turns  out  that  for  1/  =  1  the  leading  order 
correction  is  0(«:/n(K)),  [5],  [8].) 

8.3  The  transonic  layer 

Here  we  give  the  results  in  the  transonic  layer,  near  the  end  of  the  reac¬ 
tion  zone.  Note  that  we  make  an  expansion  where  the  values  of  U  and  A 
are  assumed  to  be  close  to  the  values  determined  from  the  generalized  CJ 
conditions.  Recall  that  the  values  of  their  were  determined  to  leading  order 
according  to 


1  -  A  =  where  z.  = 

7 


27^ 


and  Un  = 


7+1 


+ 


(7 + in 


7 


(205) 
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where  D'^  is  the  general  form  perturbation  of  the  detonation  velocity,  with 
out  its  order  specified. 

What  form  does  the  expansion  of  U  take  in  this  layer?  The  question  can 
be  answered  by  looking  at  the  leading  order  of  the  outer  solution  (the  near¬ 
shock  solution)  and  evaluating  it  at  the  generalized  CJ  point.  Note  that  if 
we  do  that  formally,  we  get 


7 

7-1-1  7-1-1 


(206) 


This  simple  result  strongly  suggests  the  form  on  the  TSL  expansion  must 
take.  We  introduce  the  new  layer  coordinate 


{z^Ky^‘'s  =  (1  —  A). 

Notice  that  the  term  appears  explicitly. 

We  write  the  expansion  in  the  TSL  as 

7 


(207) 


^TSL  ^ 

n  7-H 


y-l-(z./c)'/2*'u('/"*'>(s)-l- 


(208) 


The  governing  equation  for  derived  from  the  master  equation  is 


a(u(V2./))2 


ds  (7  -I- 1) 

subject  to  the  boundary  condition 

=  0  at  s  =  0 

The  solution  to  the  above  problem  is 


u 


(1/2*/)  _ 


1 


7  +  1 

The  TSL  solution  is  given  by 


A-t' 


1 


-  1 


1-1/ 


1/2 


(209) 


(210) 


(211) 


UjSL  ^ 


7  -f  1 


7 + 1 


,1-*/ 


s-  1  - 


- 1 


1-1/ 


1/2 


K 


X/lu 


(212) 
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8.4  Summary 

In  summary,  we  have  solved  for  the  solution  through  the  reaction  zone  and 
obtained  the  leading  order  corrections  that  are  0(k)  and  have  found  that 
the  solution  is  split  into  two  regions,  a  near-shock  layer  (called  the  main 
reaction  layer  MRL)  and  a  near-sonic  region,  near  complete  reaction  (the 
transonic  layer  TSL).  Importantly,  in  order  for  the  solutions  to  match,  and 
so  that  we  obtain  a  uniformly  valid  description,  we  must  have  an  eigenvalue 
relationship  between  the  perturbation  of  the  detonation  velocity.  When  we 
carry  out  an  expansion  using  the  curvature  as  our  perturbation  parameter, 
this  is  equivalent  to  finding  a  specific  value  for  the  perturbation  = 
—  Jq  F(£)di.  If  we  make  the  simplest  choice  that  the  rate  premultiplier  fc  is  a 
constant,  then  the  integral,  defined  by  the  integrand  F  in  (202)  can  be  done 
analytically  and  obtains 


=  -r 


1  27.2 


1 


+ 


1 


(213) 


fc  (7  -1- 1)2  [2(1  -  I/)  '  3  -  2i/  ■  2(2  -  I/) 

This  last  formula,  combined  with  (194),  is  the  £)„  —  k  relation  for  the  model 
explosive. 

We  now  give  a  final  set  of  (dimensional)  formulas  that  sununarize  the 
basic  result.  For  an  explosive  material  modeled  by  a  polytropic  equation  of 
state,. 


e  =  -  QK  (214) 

7-  Ip 

and  a  rate  late  of  the  form 

r  =  fc(l  -  A)*'  for  0  <  t/  <  1,  (215) 

the  (dimensional)  —  /c  relation  is  given  to  0(k)  by  the  formula 

K. 

This  representative  result  for  the  /)„  —  «:  relation  and  the  description  of  the 
reaction  zone  structure,  that  accompanies  it,  has  become  the  basis  for  the 
engineering  Method  of  Detonation  Shock  Dynamics,  that  will  be  described 
in  a  sequel  to  these  lectures. 


Dn  =  Dcj 


Dh  27^ 
fc  (7  +  I)' 


1  2  1 
+  :: — + 


2(1 -j/)  3-2i/  2(2-1/) 
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